function [RbUp,gam_p,rkhatsU,rkhatU] = find_RbU(RbU,LU,PU,nu_p,alpha_p,siga_p,lam_p,delta_p,beta_p)

Gam = 0.01:0.01:0.99;
N = length(Gam);
sigrk_p = (1+nu_p)/(alpha_p*nu_p+1)*siga_p;
ZZ = zeros(N,1);
ZZ(1) = find_gam_p(Gam(1),RbU,LU,PU,lam_p,delta_p,sigrk_p);
ii = 0;
for j=2:N
    ZZ(j) = find_gam_p(Gam(j),RbU,LU,PU,lam_p,delta_p,sigrk_p);
    if ZZ(j)*ZZ(j-1) < 0
       ii = ii+1;
       II(ii) = j;
       break
    end
end
RbUp=[];gam_p=[];rkhatsU=[];rkhatU=[];
if exist('II','var')
    gam_p = fzero('find_gam_p',[Gam(II(1)-1),Gam(II(1))],optimset('Display','on','TolFun',.1e-13,'TolX',.1e-20),RbU,LU,PU,lam_p,delta_p,sigrk_p);
    [~,rkhatsU,rkhatU] = find_gam_p(gam_p,RbU,LU,PU,lam_p,delta_p,sigrk_p);
    pi = 3.141592653589793;
    murkU = rkhatU;
    Erkd = exp(murkU+sigrk_p^2/2)*(1/2)*(1+erf((rkhatsU-murkU-sigrk_p^2)/(sigrk_p*sqrt(2))));
    RbUp = (1/beta_p - (Erkd+(1-delta_p)*PU).*LU./(LU-1))./(1-PU-lam_p.*PU);
end


end